Add SetProg templates - #38

This adds a template whose node functions are the sets of SetProg. For a homogeneous node function of degree d, the stability edge inequality V_dst(A x) ≤ γᵈ V_src(x) is the set inclusion (A / γ) S_src ⊆ S_dst, and this template hands that inclusion to SetProg.

What it adds over the existing templates is the piecewise semi-ellipsoid, a function
that is quadratic on each cone of a fixed partition. SetProg also covers the families PCC
already has, which is how the template is checked.

using SetProg
P = polyhedron(hrep([HalfSpace(a, 6.0) for a in normals]), lib)   # the partition
template = PCC.SetProgTemplate(Ellipsoid(symmetric = true, piecewise = P))
certificate = PCC.jsr_bound(template, graph, problem; optimizer)

Supported set variables

Pass one variable to share it across nodes, or a vector with one per node to give each
node its own partition. With one per node, reindex gives each node of a lifted graph
the variable of its origin, like the polyhedral templates.

Benchmarks

The system is the one from test/templates/sum_of_squares.jl: a rotation and a shear, both
scaled by 0.69. Every rate is an upper bound on the joint spectral radius from
jsr_bound with Clarabel and rtol = 1e-4. The joint spectral radius is at least 0.69,
the spectral radius of each mode.

One node (two self-loops):

de_bruijn(1, 2) (two nodes):

• Cross-checks: each of the three families PCC already has gives the same rate as its
  existing template.
• Partitions: the octagon's vertices are on the rays ±(1, 0), ±(1, 1), ±(0, 1) and
  ±(−1, 1). Its facet cones are exactly those of planar_conic_partition(1), which is
  what makes the polytope comparison meaningful.
• Piecewise semi-ellipsoids: with 8 pieces the bound drops from 0.9758 to 0.9362.
  With 4 pieces on this system, the optimum has all pieces equal to the ellipsoid.

How it is built

• Edges: add_domination! folds the scale into the map, (map / γ) * S_src ⊆ S_dst
  with γ = scale^(1/d). A nonzero margin and a scale that is a JuMP expression both
  throw an ArgumentError, as the interface asks.
• Normalization: add_normalization! constrains S ⊆ unit, i.e. V(x) ≥ ‖x‖₂ᵈ, or
  V ≥ the gauge of P for a polytope. The unit set is written in the variable's own
  family. The edge conditions are homogeneous, so this fixes the scale and excludes
  V ≡ 0 without excluding any certificate.
• Space: everything is modeled in SetProg's primal space, where a set is represented
  by its gauge function. A constant polynomial or piecewise superset has no dual
  representation.
• Fitted functions: SetProgFunction is callable. On a piecewise set it evaluates the
  piece containing x, and the function is continuous, so a boundary point can use
  either neighbour.

Limitations

• Stability only, like the other non-quadratic templates.
• No domination_slack, the same as SumOfSquaresTemplate. So edge_slacks,
  tight_edges and refine do not work with this template yet.
• No closures declared, so order_witness only uses the template-free simulation
  map. Both rows are in the graph-layer table of What works with what.
• Partitions need exact facets. Rounding noise in the facet normals, like the 1e-16
  of a cos(π/2), can make the solver report ALMOST_OPTIMAL for a point that violates
  the edge inequalities by a lot. In a SetProg prototype with Hypatia, a 45°-rotated
  square with noisy normals "certified" 0.62, below the 0.69 lower bound; with the
  normals rounded it gives 0.9758. PCC counts ALMOST_OPTIMAL as feasible, so the
  docstring warns about it and the tests use exact partitions. The cause is in
  SetProg/Polyhedra's handling of the pieces and is left for a separate fix.